package scu.maqiang.cte;

import java.io.FileWriter;
import java.io.IOException;
import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3H201;
import scu.maqiang.fes.FES3H202;
import scu.maqiang.fes.FES3H20CTE;
import scu.maqiang.mesh.Mesh3H20;
import scu.maqiang.mesh.Mesh3H8;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.mesh.Tecplot;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class BrickResultH20_DC2 {

	public static void main(String[] args) {
		double kxx = 1.0, density = 1.0, specificHeat = 1.0; 
	    double E = 2.0, Nu = 0.3, alpha = 0.05;
	    double beta = E / (1 - 2 * Nu)  * alpha;
	    
	    double theta = 0.5, T = 1.5e1, dt = T / 1500;
	    double omega = 0.25;

	    Mesh3H8 mesh3D0 = new Mesh3H8().cube(4, 4, 8);
	    mesh3D0.scale(0.5, 0.5, 1.0);

	    Mesh3H20 mesh3D = new Mesh3H20(mesh3D0);
	    
	    FES3H201 ffs = new FES3H201(mesh3D);
	    FES3H202 fsV = new FES3H202(mesh3D);
	    
	    FES3H20CTE fs = new FES3H20CTE(mesh3D, ffs, fsV);
	    
	    SRMatrix L = new SRMatrix(fs.GetNdof());
	    SRMatrix R = new SRMatrix(fs.GetNdof());
	    SRMatrix ElasticMass = new SRMatrix(fsV.GetNdof());
	    SRMatrix HTEMatrix = new SRMatrix(fsV.GetNdof(), ffs.GetNdof());
	    
	    double[] x = new double[fs.GetNdof()];
	    double[] RHS = x.clone();
	    
	    double[] temp = new double[ffs.GetNdof()];
	    double[][] dispUV = new double[3][ffs.GetNdof()];
	    
	    double[] totalDisp = new double[fs.GetNdof() / 4 * 3];
	    double[] totalDisp0 = totalDisp.clone();
	    double[] vel = totalDisp.clone();
	    double[] vel0 = totalDisp.clone();
	    double[] acc = totalDisp.clone();
	    double[] acc0 = totalDisp.clone();
	    
	    double[] ex = temp.clone();
	    double[] ey = temp.clone();
	    double[] exy = temp.clone();
	    double[] ez = temp.clone();
	    double[] tx = temp.clone();
	    double[] ty = temp.clone();
	    double[] txy = temp.clone();
	    
	    double MM = E * (1 - Nu) / ((1 + Nu) * (1 - 2 * Nu));
	    double TRef = 8.0 * 35.0 / 13.0;
	    double AA = MM + TRef * beta * beta / (density * specificHeat);
	    System.out.println("A = " + AA);
	    System.out.println("M = " + MM);
	    
	    ScalarFunc initTemp = (xy, label, param) -> xy[2] > 0.999999? 1.0:0.0;
	    double[] Itemp = ffs.valueFromFunc(initTemp, null);
	    fsV.assembleMass(new double[] {density}, BVPType.CONSIST_MASS, ElasticMass);
	    fs.assembleHeatToElastic(beta, BVPType.COMMON, HTEMatrix);
	    HTEMatrix.mul(Itemp, totalDisp);
	    NewIterSSolver solver = new NewIterSSolver(ElasticMass);
	    solver.PCGSSOR(totalDisp, acc, 1.5, 1);
	    Arrays.setAll(acc0, i -> acc[i]);
	    double[] coef = {density, specificHeat, kxx, TRef, theta, E, Nu, beta, dt, omega};
	    fs.assembleDierctCouplingMatrix_Type2(coef, BVPType.COMMON, BVPType.CONSIST_MASS,
	    		                              BVPType.COMMON, BVPType.CONSIST_MASS, BVPType.COMMON, L, R);
	    //L.sort();
	    //System.out.println(L);

	    fs.applyBCHeat_MBN(L, 2);
	    
	    fs.applyBCDisp_MBN(L, Direct.X, 3, 4, 5, 6);
	    fs.applyBCDisp_MBN(L, Direct.Y, 3, 4, 5, 6);
	    fs.applyBCDisp_MBN(L, Direct.All, 1);
	    //System.out.println(L);
	    solver.setMatrix(L);
	    double time = 0.0;
	    double[] RHS1 = new double[fs.GetNdof()];
	    
	    int lineNum = 100;
	    double dis = 1.0 / lineNum;
	    double[] lineX = new double[lineNum + 1];
	    double[] lineY = lineX.clone();
	    
	    Arrays.setAll(lineX, i -> i * dis);
	    String dir = "BrickResultH20\\";
	    String XYLineFile = dir + "Line150_2.dat";
	    Tecplot.lineXYTime(XYLineFile, lineX, lineY, 0.0);
	    
	    String heatFile = dir + "Heat150_2.dat";
	    mesh3D.toTecplot(heatFile, 0.0, x);
	    
	    String deformFile =  dir + "Deform150_2.dat";
	    mesh3D.toTecplotDeform(deformFile, 0.0, dispUV, 0.0);
	    
	    String expansionFile = dir + "expansion150_2.dat";
	    String temp00File = dir + "temp00_150_2.dat";
	    FileWriter fw1 = null;
	    FileWriter fw2 = null;
	    try {
			fw1 = new FileWriter(expansionFile);
			fw2 = new FileWriter(temp00File);
			fw1.write("0.0\t0.0\n");
			fw2.write("0.0\t0.0\n");
		} catch (IOException e) {
			// TODO Auto-generated catch block
			e.printStackTrace();
		}
	    while(time < T) {
	    	time += dt;
	    	System.out.println("time: " + time);
	    	fs.applyBCHeat_MBN(RHS, 1.0, 2);
	    	
	    	fs.applyBCDisp_MBN(RHS, Direct.X, 0.0, 3, 4, 5, 6);
		    fs.applyBCDisp_MBN(RHS, Direct.Y, 0.0, 3, 4, 5, 6);
		    fs.applyBCDisp_MBN(RHS, Direct.Z, 0.0, 1);
		    
		    solver.PCGSSSOR(RHS, x, 1.5, 1);
		    fs.extractUVH(x, temp, dispUV);
		    mesh3D.toTecplot(heatFile, time, temp);
		    System.out.println("T min: " + MVO.min(temp)  + "\t max: " + MVO.max(temp));
		    System.out.println("U min: " + MVO.min(dispUV[0]) + "\t max: " + MVO.max(dispUV[0]));
		    System.out.println("V min: " + MVO.min(dispUV[1]) + "\t max: " + MVO.max(dispUV[1]));
		    System.out.println("V min: " + MVO.min(dispUV[2]) + "\t max: " + MVO.max(dispUV[2]));
		    
		    for(int i = 0, n = ffs.GetNdof(); i < n; i++) {
		    	totalDisp[3 * i] = dispUV[0][i];
		    	totalDisp[3 * i + 1] = dispUV[1][i];
		    	totalDisp[3 * i + 2] = dispUV[2][i];
		    }
		    
		    Arrays.setAll(acc, i -> 1/(omega*dt*dt) * (totalDisp[i] - totalDisp0[i]) -1/(omega * dt) * vel0[i] - (0.5 /omega - 1) * acc0[i]);
		    Arrays.setAll(vel, i -> vel0[i] + ((1-0.5) * acc0[i] + 0.5 * acc[i]) * dt);
		    
	        Arrays.setAll(acc0, i -> 1/(omega*dt) * vel[i] + (0.5/omega -1) * acc[i]);
	        ElasticMass.mul(acc0, vel0);
	        for(int i = 0, n = fsV.GetNdof(); i < n; i++) {
	        	RHS1[i] = vel0[i];
	        }
	        R.mul(x, RHS);
	        Arrays.setAll(RHS, i -> RHS[i] + RHS1[i]);
	        Arrays.setAll(acc0, i -> acc[i]);
	        Arrays.setAll(vel0, i -> vel[i]);
	        Arrays.setAll(totalDisp0, i -> totalDisp[i]);

		    mesh3D.toTecplot(heatFile, time, temp);
		    mesh3D.toTecplotDeform(deformFile, time, dispUV, (1-Nu)/(alpha * (1 + Nu)));
		    Arrays.setAll(lineY, i -> ffs.getValue(temp, new double[] {0.25, 0.25, lineX[i]}));
	        Tecplot.lineXYTime(XYLineFile, lineX, lineY, time);
	        try {
//	        	fw1.write(time + "\t " + ffs.getValue(dispUV[2], new double[] {0.25, 0.25, 1.0}) * (1-Nu)/(alpha * (1 + Nu))  + "\n");
	        	fw1.write(time + "\t " + MVO.max(dispUV[2]) * (1-Nu)/(alpha * (1 + Nu))  + "\n");
				fw2.write(time + "\t " + ffs.getValue(temp, new double[] {0.0, 0.0, 0.0}) + "\n");
			} catch (IOException e) {
				// TODO Auto-generated catch block
				e.printStackTrace();
			}
	    }
	    
	    try {
	    	fw1.flush();
			fw1.close();
			fw2.flush();
			fw2.close();
		} catch (IOException e) {
			e.printStackTrace();
		}
	}

}
